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Abstract 

We demonstrate an application of the spectral method as a numerical approximation 
for solving Hyperbolic PDEs. In this method a finite basis is used for approximating the 
solutions. In particular, we demonstrate a set of such solutions for cases which would 
be otherwise almost impossible to solve by the more routine methods such as the Finite 
Difference Method. Eigenvalue problems are included in the class of PDEs that are 
solvable by this method. Although any complete orthonormal basis can be used, we 
discuss two particularly interesting bases: the Fourier basis and the quantum oscillator 
eigenfunction basis. We compare and discuss the relative advantages of each of these two 
bases. 
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1 Introduction 

Partial differential equations are ubiquitous in science and industry, since the dynamical laws 

governing all physical phenomena can be usually approximated by a set of partial differential 

equations. These include diverse phenomena such as the dynamics of fluids, gravitational 

fields, electromagnetic fields, etc. In particular, we are more interested in the applications in 

quantum cosmology and quantum mechanics. 
*Email: pedram@sbu.ac.ir 



1 



In some particular physical applications the problem is so simplified that the resulting 
PDE is simple enough to be solved analytically, for example by the method of separation of 
variables. However, when a more realistic modeling of the problem is required, the resulting 
partial differential equation might become so complicated that an analytical solution might 
not be feasible any more. In such cases, we have to resort to an approximate technique. 
These approximate techniques could be either analytic or numeric. An obvious example of 
an approximate analytic technique would be the usual perturbation method, in which the 
problem is divided into two segments. The main segment is supposed to be exactly solvable. 
The second segment is supposed to modify the solution obtained in the main segment only 
very slightly. This modification can be obtained analytically for any desired degree of accuracy. 
On the other hand, the numeric solutions could be either perturbative in nature or completely 
numeric. In the first explained above, the main part is solved analytically. However, 

the perturbation part is solved numerically. Examples of the completely numerical methods 
range from the simple Finite Difference Methods (FDM) [1] to the more sophisticated Multigrid 
Method, Collocation Method [2, 3], and Finite Element Method (FEM) [4, 5]. In these methods 
the configuration space is discretized, and the value of the solution at each grid point is 
determined by the values of its neighboring points. Therefore in these methods the smoothness 
of the solution is an important condition to get a reasonable approximate solution. However, 
sometimes the solutions are not smooth. Moreover, we have encountered problems which seem 
to posses intrinsic instabilities that we were not able to overcome, even when we implemented 
the Courant stability condition [6] using FDM, or the implicit FDM, or their combination. 
The simplest problem we have encountered that seems to posses both of the aforementioned 
problems is the following hyperbolic PDE. This problem appeared in applications of Robertson- 
Walker quantum cosmology with zero curvature and is a particular case of the so called the 
Wheeler-DeWitt (WD) Eq. [7], 




(1) 
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Figure 1: Left, plot of the \Reip(u, v)\ 2 of the solution to the Wheeler-DeWitt Eq. (1) with 130 
basis states used in the expansion of the solution. Right, \Keif)(u, v)\ 2 of this solution along 
the classical circular path [18] (i.e. along its crest .) 

This equation is exactly solvable [8]. By choosing an appropriate set of initial conditions, 
the absolute value squared of the solution (sometimes called the wave packet) has a smooth 
behavior and coincides well with the classical solution [8]. However, if we consider the real and 
imaginary parts of the solution separately, we can see that each part has pervasive oscillations 
almost everywhere, in particular along the crest of the \ip(u,v)\ 2 (Fig. 1). From the Figure 
it is obvious that the usage of explicit or implicit FDM would fail in this type of situations. 
Moreover for smaller classical path radii, although there will be less oscillations, the solutions 
that we attempted to find using FDM showed divergent behavior at large distances. 

We are interested in solving problems which are generalizations of the one mentioned above. 
This gives us motivation to use a different numerical method to solve this type of problems. 
This method, which was first introduced by Galerkin, consists of first choosing a complete 
orthonormal set of eigenstates of a, preferably relevant, hermitian operator to be used as 
a suitable basis for our solution. For this numerical method we obviously can not choose 
the whole set of the complete basis, as these are usually infinite. Therefore we make the 
approximation of representing the solution by only a finite superposition of the basis functions. 
By substituting this approximate solution into the differential equation, a matrix equation is 
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obtained. The expansion coefficients of these approximate solutions could be determined by 
eigenvalues and eigenfunctions of this matrix. This method has been called the Galerkin 
Method, and is a subset of the more general Spectral Method (SM) [9, 10, 11, 12]. Spectral 
methods fall into two broad categories. The interpolating, and the noninterpolating method. 
The first category, which includes the Pseudospectral and the Spectral Element Methods, 
divides the configuration space into a set of grid points. Then one demands that the differential 
equation be satisfied exactly at a set of points known as the collocation or interpolation points. 
Presumably, as the residual function is forced to vanish at an increasingly larger number of 
discrete points, it will be smaller and smaller in the gaps between the collocation points. 
The noninterpolating category includes the Lanczos tau-method and the Galerkins method, 
mentioned above. The latter is the method that we apply and, in conformity with the usual 
nomenclature, we shall simply refer to it as the Spectral Method. The interesting characteristic 
of this method is that it is completely distinct from the usual spatial integration routines, 
such as FDM, which concentrate on spatial points. In SM the concentration is on the basis 
functions and we expect the final numerical solution to be approximately independent of the 
actual basis used. That is we expect the approximate solution to converge to the exact solution 
as the number of basis elements used increases. This point has been clearly demonstrated by 
Maday et. al. [13]. Moreover in this method, the refinement of the solution is accomplished 
by choosing a larger set of basis functions, rather than choosing more grid points, as in the 
numerical integration methods. We should note that we are implicitly assuming that the 
true solution is expandable in any complete ortho normal basis such as Fourier, Laguerre[14], 
Chebyshev[15], or Legendre[16] basis. However, this requirement is usually satisfied for cases 
of physical applications. 

The paper is organized as follows: In section 2 and 3 we layout the implementation of this 
method using the quantum eigenfunctions for the simple harmonic oscillator, henceforth called 
the Oscillator basis, and the Fourier basis, respectively. In section 4 we solve a particularly 
interesting example relevant to quantum cosmology by this method using both of the afore- 
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mentioned basis functions. This problem is a generalization of the one represented in Eq. (1) 
and does not seem to have an exact solution, and we have not been able to solve this problem 
by any other numerical method that we tried. In section 5 we discuss the accuracy of each 
of the cases, and make a comparison between the two. In section 6 we discuss some general 
features of this method. 

2 The Oscillator Basis 

The general PDE equation that we want to solve is a hyperbolic one cast in the form, 

-9-2 + + ^ v - u y + ^ } ^ = °> ( 2 ) 

where f(u, v) is an arbitrary operator of u and v, but with derivatives less than two. Eq. (2) in 
general is not separable, however, any solution can be written as a superposition of the basis 
elements, 

1pm,n( U > V ) = ®m{u)Pn(v), (3) 

where each of the sets {«„} and {/3 m } is an orthonormal complete set. In this section we take 
both of them to be the Oscillator basis [17]. That is, 



where H n {x) denote the Hermite polynomials. In particular the set {ip mjn (u, v)} is a complete 
orthonormal set which can be used to span the zero sector subspace of the Hilbert space of 
the hermitian operator H, as defined in Eq. (2). These basis elements are C 2 measurable 
square integrable functions on R 2 with an inner product defined in the usual way, so that the 
orthonormality, for example, takes the form, 



J ^ m ,n{.U,v)%l) m ,, n t(u,v)dudv = 5 m ,m'5n,n' 



We can construct a general solution as follows, 



i/;(u, v) = J2 A min a m (u)f3 n (v), (6) 

m,n 

and can use the following expansion, 

f(u, V)lp(ll, V) = J2 B m,nOt m {u)l3 n {v), (7) 
m,n 

where B m n are coefficients that can be determined once f(u, v) is specified. By substituting 
Eqs. (6,7) in Eq. (2), and using the differential equation of the Hermite polynomials we obtain, 

]T [[(2m + l)wi - (2n + l)cu 2 ] A m>n + B m , n ] a m {u)(5 n {v) = 0. (8) 

m,n 

Because of the linear independence of a m (u)s and /3 n (v)s, every term in the summation must 
satisfy, 

[(2m + l)u! - (2n + l)uJ 2 }A m , n + B m>n = 0. (9) 
It only remains to determine the matrix B. By using Eqs. (6,7) we have, 

B m,na m (u)(3 n (v) = Y An,n/(«, v)a m (u) (3 n {v) . (10) 
m,ra m,n 

By multiply both sides of the above equation by a m i(u)/3 n i(v) and integrating over the full 
range of variables u and v, and using ortho normality of the basis functions, one finds, 

/ / a m (u)(3 n (v)f(u, v)a m ,{u)(3 nl (v)dudv ) A m , , a , = C ™ 

,n,m' ,n' A m ' jU ' . 

n, .„ J - ooJ - 00 / m ',n> 

(11) 

Therefore we can rewrite Eq. (9) as, 

[(2m + 1)^! - (2n + l)w 2 ] An,n + H C m,n,m',n' An'.n' = 0. (12) 

m',n' 

It is obvious that the presence of the operator /(«, u) in Eq. (2), leads to nonzero coefficients 
C m ,n,m',n' in Eq. (12), which in principle could couple all of the matrix elements of A. For 
the usual choices of f(u,v), e.g. the choice presented in [8]: / = \k(u 2 — v 2 ) 1 / 3 , the problem, 
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to the best of our knowledge, is not analytically solvable. Therefore we have to resort to a 
numerical solution. In general the number of basis elements are at least countably infinite. 
The aforementioned coupling of terms in the main matrix Eq. (12) forces us to make the 
approximation of using a finite basis. It is easy to see that the more basis functions we 
include, the closer our solution will be to the exact one. We select the first N basis functions 
in each direction, that is m and n run from 1 to N. Then we replace the square matrix A with 
a column vector A' with N 2 elements, so that any element of A corresponds to one element of 
A 1 . With this replacement, Eq(12) can be written as, 

DA' = 0, (13) 

where D is a square matrix with N 2 x N 2 elements which can be obtained from Eq. (12). 
Looked upon as an eigenvalue equation, i.e. DA' a = aA' a , the matrix D has N 2 eigenvectors. 
However, for constructing the acceptable wavef unctions, i.e. the ones satisfying the WD Eq. 
(2), we only require eigenvectors which span the null space of the matrix D. That is, due to Eq. 
(12) we will have exactly iV null eigenvectors which will be linear combination of our original 
eigenfunctions introduced in Eq. (3). After finding these N eigenvectors A' 1 (i = 1, 2, 3, A), 
we can find the corresponding elements of the matrix A, A l mn . Therefore the wavefunction 
can be expanded as, 

ij(u,v) = J2^¥(u,v), (14) 

i 

where ip l (u,v) are the appropriate eigenfunctions for the problem {i.e. Eq. (2)), 

^( U , V) = J2 ^ m ,n^m{u)l3 n {v). (15) 
m,n 

In Eq. (14) A*s are arbitrary complex constants to be determined by the initial conditions. 

3 The Fourier Basis 

As mentioned before, any complete orthonormal set can be used for this SM. In this section 
we use the Fourier series basis as a second example. That is, since we need to choose a finite 
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subspace of a countably infinite basis, we restrict ourselves to a finite square region of sides 
2L. This means that we can expand the solution as, 

2 



. . . ^ ^ . /miTu\ fnirv\ 

1 f J (u, v)= 2^ 2^ A m ^ gi —— 9j[—j-) 

i,j=lm,n \ L / \ L J 



(16) 



where, 



* ( H f H ) = v / S sin ( I f M )' 

92 ( H f H ) = V / S C0S ( ! T)- 



and R % 



1, m ^ 

2, m = 



That is in the Fourier basis we assume periodic boundary condition. By referring to the WD 
Eq. (2), we realize that in the Fourier basis it is appropriate to introduce /' as, 



f'[u, v) = f(u, v) + ufu 2 - u? 2 v 2 , 
and, in analogy with Eq. (7), we can make the following expansion, 

i,j m,n \ L J \ L 

By following steps analogous to those of Eqs. (8-11) we obtain, 



(17) 



(18) 



mix 



rnr 

T 



A m: n,i,j + B mni j 0, 



(19) 



where 



B' 

JD m,n,i,j 



E 



m',n' 



L ( miru\ frnrv \ f.. . / m'nu\ ( n'7rv\ 

l 9i \ T~ ) 9j ) ^' v ' 9i ' \ _ L~ I 9i ' { ~T~ j 



A m ' t n',i' ,j' 



E ^m,n,i,j,m',n',i',j' An' ,ri 
m',n',i',j' 

Therefore we can rewrite Eq. (19) as 



(20) 



mnf - (mr) 2 ] A m>nAd + ^ C' mMm , ^ ^ jf An',n',*'j' = 0. 



(21) 



In this case, we select 4iV 2 basis functions. Using reasoning analogous to the previous case, 
for example by defining the column vector A' out of the matrix A, we transform Eq. (21) to 



D' A' — 0. 



(22) 
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Where, as before, D' is a square matrix now with (2N) 2 x (2iV) 2 elements which can be easily 
obtained from Eq. (21). After finding the 4N 2 eigenvectors of D', we select the 2N ones with 
zero eigenvalue, i.e. A' k (k = 1, 2, 3, 27V). We can then find the corresponding elements of 
matrix A, A k mni y Therefore, the wavefunction can be expanded as 

iP(u,v) = '£\ k iJ> k (u,v), (23) 
k 

where, as before, t/j k (u,v) are the appropriate eigenfunctions for the problem (i.e. Eq. (2)), 

m,n,i,j 

Here A fe s in Eq. (23) are again arbitrary complex constants to be determined by the initial 
conditions. 

Now we apply this method to one of the examples stated above which happens to be 
relevant in quantum cosmology, and was our original motivation for using this method. 

4 Application of the Spectral Method to a Specific Ex- 
ample 

For a specific example, we consider a hyperbolic PDE which happens to be the Wheeler-DeWitt 
equation for the Robertson- Walker quantum cosmology with non-zero curvature, 

~& 2+ & 2+ ^ - ^ + 1^ " V2)1/3 } ^ V) = °- (25) 

As mentioned before, the case k = is exactly solvable [8] and has a closed form solution in 
the Oscillator basis. We shall state these solutions here for illustrative purposes, especially for 
motivating the choice of initial conditions and comparison with the non-trivial cases, i.e k ^ 0, 
which shall be solved by this method next. In this case f(u,v) = in Eq. (2), so the matrix 
B in Eq. (9) is also zero due to Eq. (7). Therefore, Eq. (9) reduces to, 

[(2m + 1 Vx - (2n + l)w 2 ] A„ hn = 0. (26) 
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This means A mn = for, 



/ 2m + 1 , nn . 

U2 * ^TT Ul - (27) 



That is we have nontrivial solutions only for u\ and u>2 being a rational multiple of each other. 
Choosing uj\ = uj 2 for simplicity, we have n = m and the expansion of the wavefunction (Eq. 
(6)) reduces to, 

^(u, v) =J2 A m a m (u)f3 m (v), (28) 

m 

where A rn = A m ^ m are complex expansion coefficients that can be determined by applying ini- 
tial conditions (i.e. specifying ip(u,0) and ^\ v=0 ). In references [8, 18] the authors considered 
the following initial condition on the wave function, 

0) = (e-M 2 /2 + e -(x+x)V2) . (29) 

We can decompose any initial wavefunction in the Oscillator basis. For the choice presented 
in Eq. (29) we have, 

VKw,0) = J2 c n a n(u), where c n = e " lx|2 / - — = , (30) 
n V2 n n! 

and x is an arbitrary complex number. The prime on the summation denotes the restriction 
that n is even. These coefficients are same as those of the coherent states of a one dimensional 
simple harmonic oscillator. With this choice of coefficients we expect that a classical-quantum 
correspondence should be manifest. The canonical choice for initial slope is [8], 



dv 



\- c n a n (u)H' n (®) 

(_!)(n/2) n | ' \ 6L ) 



v=0 n (V2)T 

where H n s are the Hermite polynomials and prime denotes differentiate respect to v. By 
choosing % to be a real, the classical paths corresponding to these solutions can be shown 
to be circles with radii x- We have found that 35 basis functions are sufficient for finding 
the solution to an accuracy of about 10~ 8 , when the classical radius of the wave packet (x) 
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Figure 2: Oscillator Basis: Left, the absolute value squared of the wave packet \if)(u,v)\ 2 for 
X = 4 and N = 35, k = 0. Right, the contour plot of the same figure with the classical path 
superimposed as the thick solid line. 

is less than 4 (Fig. 2). As can be seen in the figure, and also for all the cases presented in 
[8], the classical-quantum correspondence is manifest. Note that the wave packet has compact 
support and in the oscillator basis, the truncation of the basis functions automatically restricts 
the solution to have this property. We only need to choose the configuration space domain to 
be large enough. 

At this point it seems to us that a brief mention of the relevant dynamical equations for the 
classical cosmology might be helpful, at least for completeness, especially in light of the fact 
that "classical-quantum correspondence" that we keep mentioning in this paper is particularly 
important in physics. Moreover, the non-linearity and the moving singular behavior of these 
equations would become apparent. These equations are (see, for example, [19]), 

3k u 

u + u + T (u 2 -v 2 ) 2 /* = °' (32) 

3 Aj u 

" + v + Y (u 2 -v 2 ) 2 ^ =^ (33) 

u 2 + u 2 - v 2 - v 2 + ^k(u 2 - v 2 ) 1 / 3 = 0. (34) 
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Here the u and v variables are functions of time. Equations (32) and (33) are the dynamical 
equations and, Eq. (34) is the zero energy constraint, from which the Wheeler-deWitt Eq. (25) 
arises. An appropriate initial conditions is the following, 

«(0) = -X, u(0)=0, is(0)=0, i)(0) = i)o, (35) 

where x can be treated as a free parameter and v is adjusted so that Eq. (34) is satisfied. 
With this choice of initial conditions the classical path would be, as mentioned before, exactly 
a circle in the k = depicted in Fig. 2. In the general case of ou\ ^ u 2 , these solutions 

would give Lissajous figures. 

For the case k ^ 0, the problem is not exactly solvable in quantum cosmology and we 
will use the SM to get an approximate solution. We have to mention that the corresponding 
equations for classical cosmology, Eqs. (32)-(34), are a set of non-linear, coupled ODEs with 
moving singularities which are not exactly solvable either. However, a general method for 
solving them has been presented in [20], and this is the method we shall use. Also a detailed 
explanation of the physical setting of the problem in the classical domain has been presented 
in [19]. On one hand, we need to choose a set of appropriate initial conditions for both the 
classical and quantum cases which would make their correspondence manifest when they are 
superimposed, as in Fig. 2 which was for the k = case. However, we should note that the 
values of v in the classical case depends on k. On the other hand, an appropriate choice for 
the initial conditions should be such that we could easily compare our results with the k = 
case. Therefore, here we choose the same initial conditions for the quantum cosmology Eqs. 
(30,31), and classical cosmology Eq. (35), as for the k = case. An alternative would be to 
choose appropriate initial canonical slope for the quantum cosmology case [21]. Having set the 
initial conditions, we proceed to solve the quantum cosmology problem for the case k ^ with 
the Oscillator basis. 

Both the classical and quantum solutions for the case k = +1 are shown in Fig. 3. We can 
see that the general behavior of this case is similar to the k = case. In particular, again we 
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Figure 3: Oscillator Basis: Left, the absolute value squared of the wave packet \if)(u,v)\ 2 for 
X = 4 and N = 35, k = +1. Right, the contour plot of the same figure with the classical path 
superimposed as the thick solid line. 

have a very good classical-quantum correspondence. However, although the extreme points of 
the solution do not change as compared to the k = case, the whole pattern is a little wider. 
Moreover, the solution (l^l 2 ) is not as smooth as k = case. 

Both the classical and quantum solutions for the case k = —1 are shown in Fig. 4. We can 
see that the general behavior of this case is also similar to the k = case. In particular, again 
we have a very good classical-quantum correspondence. However, although the extreme points 
of the solution do not change as compared to the k = case, the whole pattern is a little 
narrower. Moreover, the solution (\ip\ 2 ) is not as smooth as k = case. More importantly, 
when k = — 1 we encounter a new characteristic of the solution. Note that the u 2 — v 2 term 
in Eq. (25) is usually dominant and causes stability of the solutions, for all k. However, the 
\k{u 2 — -y 2 )^/ 3 ) term in that equation could cause instability near u = ±v lines, only in the 
case of k < 0. However this is precisely where the dominant term vanishes. Therefore we do 
expect numerical instabilities along these two lines for this case. These instabilities can be 
seen in Fig. 4. However, due to numerical approximations made, the instabilities seem to be 
more pronounced along the line u = v in our solution. We have not been able to the pinpoint 
the exact source of this numerical asymmetry in the instability. 
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Figure 4: Oscillator Basis: Left, the absolute value squared of the wave packet \if)(u,v)\ 2 for 
X = 4 and N = 35, k = — 1. Right, the contour plot of the same figure with the classical path 
superimposed as the thick solid line. 



It would be an interesting comparison to solve exactly the same problem in the Fourier 
basis. By using the procedure mentioned in section 3, and again choosing exactly 35 basis 
functions for ease of comparison, we easily find the solutions for k = 0, k = 1 and k = — 1, 
which are shown in Figs. (5-7), respectively. 

5 Comparison of the Oscillator and Fourier Bases 

Here, we compare the results of the Spectral Method using the two finite bases. To this end, 
we first discuss the errors of the solutions. In the case k = we have the exact solutions as 
an infinite series. For small radii, e.g. x — 4, the difference between the exact solution and 
a series solution with 120 oscillator terms is absolutely negligible. Therefore, we can safely 
substitute this finite series for the exact one. For computing errors, we divide the 2D base 
domain into M 2 grid points. Then, we average the square of the absolute value of the difference 
of the exact solution with that obtained by the 35 Fourier or oscillator basis functions on the 
grid points: 

, 2 ES^35(^V^1 2 o(^)| 2 .... 

°k=0 — v^Af I / / • -M2 ' v ' 
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Figure 5: Fourier Basis: Left, the absolute value squared of the wave packet \ip(u, v)\ 2 for 
X = 4 and N = 35, k = 0. Right, the contour plot of the same figure with the classical path 
superimposed as the thick solid line. 




Figure 6: Fourier Basis: Left, the absolute value squared of the wave packet \ip(u, v )| 2 for 
X = 4 and N — 35, k — 1. Right, the contour plot of the same figure with the classical path 
superimposed as the thick solid line. 
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Figure 7: Fourier Basis: Left, the absolute value squared of the wave packet \if)(u, v)\ 2 for 
X = 4 and N = 35, k = —1. Right, the contour plot of the same figure with the classical path 
superimposed as the thick solid line. 

In the case k ^ the problem is not exactly solvable, so for comparison purposes we define 
a measure for the error as the average square of the absolute value of the difference of the 
solutions with, for example, 35 and 30 basis functions. That is, 





r2 




- tp3o{hj)\ 2 






i,j)\ 2 




$ Fourier 


^oscillator 


k = 


= 


6.70875 x 10~ 2 


4.16459 x 10~ 7 


k = 


+ 1 


4.01423 x 10~ 3 


3.84836 x 10~ 2 


k = 


-1 


2.83746 x 10~ 3 


3.64503 x 10~ 2 



(37) 



Table I. Errors for the Oscillator and Fourier basis 

From Table I we can conclude that the oscillator basis is more appropriate than Fourier basis 
for the case k = 0. It seems that this is only due to the fact that the Oscillator basis is the 
exact solution of the problem in this case. However, the Fourier basis is more appropriate for 
the case k ^ 0. The main reason for this seems to us to be the fact that the solutions have 
compact support, so we have an extra parameter that we can adjust (the size of spatial domain 
(2L)) and this yields better results. In fact one can use this freedom to set up an optimization 
procedure [22]. 
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6 Discussion 



Here we have exhibited the implementation of the SM (Galerkin) for solving hyperbolic PDEs. 
We use finite basis of Oscillator or Fourier eigenfunctions, for example, and show that in some 
cases where the popular numerical methods such as FDM or FEM fail, this method gives 
reasonable results very easily. The requirement that the solution should be expandable in 
a complete orthonormal basis is crucial. However, certain bases might be more appropriate 
for a given problem. For example, if the wave packets of previous section was not damped 
strongly in u and v directions, choosing oscillator basis might not have been as appropriate. We 
have found, much to our surprise, that the main source of error is the numerical integrations, 
as compared to varying the number of basis elements. This arise due to the fact that in 
order to find the coefficients C(m,n,m' ,n'), we need to calculate N 4 two fold integrations, 
where iV denotes the number of basis elements. This is the most time consuming part of 
the procedure. Using programs with refined integration routines such as Mathematica would 
have consumed too much time even at their default levels. So we used the simple trapezoid 
integration technique using Fortran, and this severely limited our accuracy when we refrained 
from refining our mesh excessively to avoid spending too much time. However, calculations of 
C(m,n,m' ,n')s can be parallelized because these coefficients are independent. In some cases 
we can calculate some of the integrals analytically (e.g. the values of C' m n i j m / in Eq. 
(20) when k — 0), which decreases the computation time. The memory consuming part of the 
procedure is solving the matrix equation (Eq. (13)) and this increases proportional to N A . 
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